2026-03-15
=====================================================
머릿말
Note. 참고 URL
https://github.com/seejokim1/NA_with_python
21세기는 기술 혁신과 지식의 융합이 핵심 동력으로 작용하는 시대다. 특히 4차 산업혁명과 인공지능 기술의 발전은 전통적인 공학 분야에 새로운 도전과 기회를 제공하고 있다. 이러한 환경에서 수치해석은 공학 문제를 해결하는 핵심 도구로, 다양한 데이터 기반 분석과 최적화 문제 해결에 활용되고 있다.
현재 수치해석의 중요성은 더욱 강조되고 있으며, 특히 자율주행 자동차의 수치해석에서 가장 중요한 것은 정확도와 실시간 처리이다. 자율주행 차량은 다양한 센서 데이터를 실시간으로 처리하며, 도로 상황을 정확하게 예측하고 반응해야 하기 때문이다. 이러한 기술은 수치해석을 통해 더욱 정교해지고, 안전한 자율주행을 가능하게 한다.
또한 AI와 수치해석은 깊은 관계를 맺고 있으며, AI의 발전이 수치해석 방법론에 많은 영향을 미쳤고, 반대로 수치해석 기술이 AI의 성능을 향상시키는 데 중요한 역할을 해왔다. 수치해석은 실세계의 문제를 수학적으로 모델링하고 해결하기 위한 다양한 알고리즘과 방법을 제공하며, AI는 이러한 수학적 모델을 데이터 기반으로 학습하고 최적화하는 능력을 갖추고 있다. 이 둘의 결합은 현대 과학과 공학의 중요한 연구 분야로 자리 잡았다.
예를 들어, 2024년 노벨 물리학상은 머신러닝과 AI 기술이 물리학 문제 해결에 중요한 기여를 한 연구자에게 수여되었다. 이는 머신러닝이 복잡한 물리적 시스템을 시뮬레이션하거나 예측하는 데 사용되고 있음을 보여준다. 또한 알파고에 이어 알파폴드와 같은 AI 기술이 노벨 화학상을 수상하며, 수치해석 기법을 바탕으로 AI 모델이 실세계 문제를 해결하는 데 어떻게 활용될 수 있는지를 보여주었다.
이 책은 파이썬 프로그래밍 언어를 활용하여 수치해석의 기초부터 심화 내용까지 단계적으로 학습할 수 있도록 구성되었다. 첫 장에서는 파이썬의 기본 문법과 데이터 처리에 필수적인 라이브러리(예: NumPy, Matplotlib)를 다루며, 프로그래밍과 수치해석의 연결 고리를 제공한다. 이어지는 장에서는 비선형 방정식 근사법, 선형 연립방정식 해법, 수치미분 및 수치적분, 보간법, 수치 회귀 등 전통적인 수치해석의 핵심 주제들을 체계적으로 소개한다.
특히 후반부에서는 상미분 방정식과 편미분 방정식의 수치적 해법을 다루며, 유한차분법(Finite Difference Method)과 유한요소법(Finite Element Method)의 응용을 통해 실질적인 문제 해결 방법을 제시한다. 이와 더불어 Python을 활용하여 다양한 수치해석 문제를 직접 해결해보는 연습문제를 포함해 독자가 이론과 실습을 병행하며 학습할 수 있도록 하였다.
김시조
유한요소법(Finite Element Method, FEM)은 복잡한 구조와 물리 시스템을 수치적으로 해석하기 위한 방법이다.
연속적인 물리 문제를 작은 요소(element)로 분할하고, 각 요소에서 근사해를 구성하여 전체 해를 조합한다.
FEM은 다음과 같은 분야에서 널리 사용된다.
구조 해석
열전달 해석
전자기장 해석
유체 해석
다물리(multiphysics) 문제
—
CAE 환경에서 FEM 해석은 다음 과정을 포함한다.
문제 정의 및 수식화
연속적인 물리 현상을 지배 방정식(PDE)으로 정의한다.
\[\mathcal{L}(u) = f\]
이를 변분형(Variational Form) 또는 약형(Weak Form)으로 변환하여 수치해석이 가능하도록 한다.
—
Shape Function 및 근사 해법
각 요소 내부에서 shape function을 사용하여 해를 근사한다.
\[u(x) \approx \sum_{i=1}^{n} N_i(x) u_i\]
여기서 \(N_i(x)\)는 shape function, \(u_i\)는 절점 자유도이다.
Galerkin 방법을 사용하여 근사 방정식을 구성한다.
—
전역 행렬 구성
요소 행렬을 조립(assembly)하여 전역 시스템을 구성한다.
\[[K]\mathbf{u} = \mathbf{f}\]
구조 문제에서는 강성 행렬, 열전달 문제에서는 전도 행렬이 형성된다.
—
경계 조건 적용
Dirichlet 또는 Neumann 조건을 전역 행렬에 반영한다.
—
해 계산 및 후처리
선형 또는 비선형 시스템을 풀어 해를 계산하고, 응력, 변형률, 온도 분포 등을 후처리한다.
—
복잡한 형상에 적용 가능
비정형(unstructured) 격자 사용 가능
다물리 문제 확장 용이
고차 요소 사용 가능
CAE 핵심 해석 기법
—
항공우주 구조 해석
자동차 충돌 해석
반도체 열해석
에너지 시스템 설계
바이오메디컬 응용
HPC와 AI 기술의 발전으로 대규모 FEM 해석의 정확도와 효율성이 지속적으로 향상되고 있다.
FEM의 기본 아이디어는 연속적인 경계값 문제(BVP)를 이산(discrete) 시스템 방정식으로 변환하는 것이다.
연속체를 여러 개의 요소(element)로 분할하고, 각 요소의 강성 관계를 구성한 후 전역 시스템으로 조립(assembly)한다.
—
각 절점에서 변위와 힘을 각각
\[u_i, \quad F_i\]
로 정의한다.
—
첫 번째 요소에 대해
\[\begin{bmatrix} k_{11}^{(e=1)} & -k_{12}^{(e=1)} \\ -k_{21}^{(e=1)} & k_{22}^{(e=1)} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \begin{bmatrix} F_1^{(e=1)} \\ F_2^{(e=1)} \end{bmatrix} \tag{12.1-1}\]
두 번째 요소에 대해
\[\begin{bmatrix} k_{11}^{(e=2)} & -k_{12}^{(e=2)} \\ -k_{21}^{(e=2)} & k_{22}^{(e=2)} \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} F_2^{(e=2)} \\ F_3^{(e=2)} \end{bmatrix} \tag{12.1-2}\]
—
요소 행렬을 전역 좌표계로 조립하면
\[[K]\mathbf{u} = \mathbf{F}\]
전역 강성 행렬은 다음과 같은 구조를 갖는다.
\[\begin{bmatrix} k_{11}^{(1)} & -k_{12}^{(1)} & 0 & 0 \\ -k_{21}^{(1)} & k_{22}^{(1)}+k_{11}^{(2)} & -k_{12}^{(2)} & 0 \\ 0 & -k_{21}^{(2)} & k_{22}^{(2)}+k_{11}^{(3)} & -k_{12}^{(3)} \\ 0 & 0 & -k_{21}^{(3)} & k_{22}^{(3)} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} = \begin{bmatrix} F_1 \\ F_2 \\ F_3 \\ F_4 \end{bmatrix} \tag{12.1-3}\]
—
내부 절점 \(j\)에 대해
\[k_{21}^{(e=j-1)} u_{j-1} + \left( k_{22}^{(e=j-1)} + k_{11}^{(e=j)} \right) u_j - k_{12}^{(e=j)} u_{j+1} = F_j\]
이는 1차원 스프링-질량 시스템의 평형 방정식과 동일한 형태이다.
—
각 절점에서 다음 중 하나를 지정할 수 있다.
변위 (Dirichlet 조건)
힘 (Neumann 조건)
적절한 경계 조건이 없으면 전역 행렬은 특이 행렬(singular matrix)이 된다.
\[\det[K] = 0\]
경계 조건을 적용하면 well-posed problem이 된다.
—
요소 행렬 구성
전역 행렬 조립
경계 조건 적용
선형 시스템 해법
FEM의 핵심은 local → global 조립 과정이다.
형상함수(Shape Function)는 요소 내부의 해를 절점 값으로 근사하기 위해 사용되는 함수이다.
일반적으로
\[u(x) \approx \sum_{i=1}^{n} N_i(x) u_i\]
여기서
\(N_i(x)\) : 형상함수
\(u_i\) : 절점 자유도
형상함수는 다음 조건을 만족한다.
\[N_i(x_j) = \begin{cases} 1 & (i=j) \\ 0 & (i \neq j) \end{cases} \tag{12.2.1-1}\]
—
기준 요소(reference element) \(\xi \in [-1,1]\)에서
\[N_1(\xi) = \frac{1-\xi}{2}, \qquad N_2(\xi) = \frac{1+\xi}{2} \tag{12.2.1-2}\]
보간식은
\[u(\xi) = N_1(\xi) u_1 + N_2(\xi) u_2 \tag{12.2.1-3}\]
또는
\[u(\xi) = u_1 \frac{1-\xi}{2} + u_2 \frac{1+\xi}{2}\]
좌표 보간 역시 동일한 형식을 따른다.
\[x(\xi) = N_1(\xi)x_1 + N_2(\xi)x_2\]
이를 등매개변수(isoparametric) 개념이라 한다.
—
노드가 3개인 경우 (2차 요소)
\[N_1(\xi) = \frac{1}{2}\xi(\xi-1)\]
\[N_2(\xi) = 1 - \xi^2\]
\[N_3(\xi) = \frac{1}{2}\xi(\xi+1)\]
2차 보간식:
\[u(\xi) = \sum_{i=1}^{3} N_i(\xi) u_i\]
—
기준 좌표 \((\xi,\eta) \in [-1,1]\times[-1,1]\)
4절점 사각형 요소에서
\[N_1(\xi,\eta) = \frac{1}{4}(1-\xi)(1-\eta)\]
\[N_2(\xi,\eta) = \frac{1}{4}(1+\xi)(1-\eta)\]
\[N_3(\xi,\eta) = \frac{1}{4}(1+\xi)(1+\eta)\]
\[N_4(\xi,\eta) = \frac{1}{4}(1-\xi)(1+\eta) \tag{12.2.1-5}\]
보간식은
\[u(\xi,\eta) = \sum_{i=1}^{4} N_i(\xi,\eta) u_i\]
—
기준 좌표 \((\xi,\eta,\zeta)\in[-1,1]^3\)
8절점 육면체 요소에서
\[N_i(\xi,\eta,\zeta) = \frac{1}{8} (1+\xi\xi_i) (1+\eta\eta_i) (1+\zeta\zeta_i)\]
여기서
\[(\xi_i,\eta_i,\zeta_i) = (\pm1,\pm1,\pm1)\]
각 절점에 대해
\[N_1 = \frac{1}{8}(1-\xi)(1-\eta)(1-\zeta)\]
\[N_2 = \frac{1}{8}(1+\xi)(1-\eta)(1-\zeta)\]
\[N_3 = \frac{1}{8}(1+\xi)(1+\eta)(1-\zeta)\]
\[N_4 = \frac{1}{8}(1-\xi)(1+\eta)(1-\zeta)\]
\[N_5 = \frac{1}{8}(1-\xi)(1-\eta)(1+\zeta)\]
\[N_6 = \frac{1}{8}(1+\xi)(1-\eta)(1+\zeta)\]
\[N_7 = \frac{1}{8}(1+\xi)(1+\eta)(1+\zeta)\]
\[N_8 = \frac{1}{8}(1-\xi)(1+\eta)(1+\zeta) \tag{12.2.1-6}\]
—
보간 조건 만족 (Kronecker delta property)
\(\displaystyle \sum_i N_i = 1\)
등매개변수 개념 적용 가능
고차 요소로 확장 가능
이 절에서는 실제 좌표계 \(x\)에서 유한요소법(FEM)을 적용하는 절차를 정리한다.
연속 문제를 약형(Weak Form)으로 변환한 후 요소 단위로 이산화하여 전역 시스템을 구성한다.
지배 방정식
다음 2차 미분 방정식을 고려한다.
\[- k \frac{d^2 u(x)}{dx^2} + c \frac{du(x)}{dx} + b u(x) = f(x), \quad x \in [0,1]\]
Dirichlet 경계조건:
\[u(0)=0, \qquad u(1)=0\]
약형(Weak Form) 유도
테스트 함수 \(v(x)\)를 곱하고 적분하면
\[\int_0^1 v(x) \left( - k u'' + c u' + b u - f \right) dx = 0\]
부분적분을 적용하면
\[\int_0^1 k u' v' + c v u' + b u v \, dx = \int_0^1 f v \, dx + \left[ k v u' \right]_0^1\]
Dirichlet 조건이 적용되면 경계항은 소멸한다.
따라서 약형은
\[\int_0^1 \left( k u' v' + c v u' + b u v \right) dx = \int_0^1 f v dx\]
형상함수 적용
각 요소 \(e=[x_1,x_2]\)에서 선형 형상함수:
\[\phi_1(x)=\frac{x_2-x}{h}, \qquad \phi_2(x)=\frac{x-x_1}{h}\]
도함수:
\[\phi_1'=-\frac{1}{h}, \qquad \phi_2'=\frac{1}{h}\]
근사해:
\[u(x) = \sum_{i=1}^{2} U_i \phi_i(x), \quad v(x)=\phi_j(x)\]
요소 강성행렬과 하중벡터
요소 방정식:
\[K^{(e)} U^{(e)} = F^{(e)}\]
강성행렬:
\[K^{(e)}_{ij} = \int_{x_1}^{x_2} \left( k \phi_i' \phi_j' + c \phi_j \phi_i' + b \phi_i \phi_j \right) dx\]
계산 결과:
\[K^{(e)} = \begin{bmatrix} \frac{k}{h} + \frac{b h}{3} - \frac{c}{2} & -\frac{k}{h} + \frac{b h}{6} + \frac{c}{2} \\[6pt] -\frac{k}{h} + \frac{b h}{6} - \frac{c}{2} & \frac{k}{h} + \frac{b h}{3} + \frac{c}{2} \end{bmatrix}\]
하중 벡터:
\[F_i^{(e)} = \int_{x_1}^{x_2} f(x)\phi_i(x) dx\]
전역 조립
전체 시스템:
\[[K]\mathbf{U}=\mathbf{F}\]
요소 행렬을 전역 행렬로 조립한 후 경계 조건을 반영한다.
핵심 절차 요약
PDE 정의
Weak Form 유도
Shape Function 적용
요소 행렬 계산
전역 조립
경계 조건 적용
선형 시스템 풀이
이 예제는 식 (12.3-10), (12.3-11)을 기반으로 1차원 FEM을 실제 좌표계에서 구현한 것이다.
지배 방정식
\[- \frac{d^2 u(x)}{dx^2} + 0 \cdot \frac{du(x)}{dx} + 0 \cdot u(x) = f(x), \qquad x \in [0,1]\]
여기서
\[f(x) = -e^x\]
—
경계 조건
Dirichlet 경계조건:
\[u(0)=1, \qquad u(1)=e\]
—
행렬 구성
요소 강성행렬은
\[K^{(e)} = \begin{bmatrix} \frac{1}{h} & -\frac{1}{h} \\ -\frac{1}{h} & \frac{1}{h} \end{bmatrix}\]
전역 조립 후
\[[K]\mathbf{U}=\mathbf{F}\]
—
경계 조건 처리
경계 조건은 강성 행렬과 하중 벡터에 직접 반영한다.
\[K[0,:]=0, \quad K[0,0]=1\]
\[K[-1,:]=0, \quad K[-1,-1]=1\]
우변 벡터는
\[F[0]=1, \qquad F[-1]=e\]
경계 외부 항은 0으로 설정한다.
—
최종 선형 시스템
\[[K]\mathbf{U}=\mathbf{F}\]
이를 풀어 절점 해를 구한다.
—
핵심 요약
Weak form 기반 요소 행렬 계산
전역 행렬 조립
Dirichlet 조건 직접 대입
선형 시스템 해법 적용
FEM에서 좌표 변환은 기준 요소(reference element)를 기반으로 계산을 표준화하고, 실제 요소의 기하학적 복잡성을 효율적으로 처리하기 위해 사용된다.
좌표 변환을 통해 복잡한 실제 형상을 단순한 기준 요소로 매핑(mapping)하여 계산한다.
기준 좌표계 변환의 필요성
복잡한 요소 형상 처리 가능
통일된 수학적 표현 사용
수치 적분(Gaussian Quadrature) 적용 용이
고차 형상 함수 확장 가능
계산 효율성 향상
기준 좌표계 (Reference Coordinate System)
기준 요소는 단순한 표준 영역에서 정의된다.
\[\xi \in [-1,1]\]
1차원: 선분 2차원: 삼각형 또는 사각형 3차원: 육면체 등
실제 좌표계 (Physical Coordinate System)
실제 요소는 물리적 영역 \(x\)에서 정의된다.
각 요소는 고유한 길이와 형상을 가진다. 실제 좌표에서 직접 적분하는 것은 복잡하므로 기준 좌표로 변환하여 계산한다.
1차원 좌표 변환
기준 좌표:
\[\xi \in [-1,1]\]
실제 좌표:
\[x \in [x_1,x_2]\]
등매개변수 변환:
\[x(\xi) = \sum_{i=1}^{2} \phi_i(\xi) x_i\]
선형 형상함수:
\[\phi_1(\xi)=\frac{1-\xi}{2}, \qquad \phi_2(\xi)=\frac{1+\xi}{2}\]
따라서
\[x(\xi)=x_1 \phi_1(\xi)+x_2 \phi_2(\xi)\]
Jacobian
\[dx = J d\xi\]
\[J = \frac{dx}{d\xi} = \frac{x_2 - x_1}{2}\]
적분 변환
\[\int_{x_1}^{x_2} f(x) dx = \int_{-1}^{1} f(x(\xi)) J d\xi\]
즉,
\[\int_{x_1}^{x_2} f(x) dx = \int_{-1}^{1} f(x(\xi)) \frac{x_2-x_1}{2} d\xi\]
핵심 정리
실제 좌표 → 기준 좌표 매핑
Jacobian을 이용한 적분 변환
Gaussian 적분은 기준 좌표에서 수행
FEM 계산의 일반성 및 효율성 확보
지배 방정식
다음 2차 미분방정식을 고려한다.
\[- a(x)\frac{d^2 u(x)}{dx^2} + b(x)\frac{du(x)}{dx} + u(x) = f(x), \qquad x \in [0,1] \tag{12.4.2-1}\]
경계 조건:
\[u(0)=0, \qquad u(1)=g(1)\]
FEM은 연속 문제를 이산 문제로 변환하여
\[KU = F\]
형태의 선형 시스템을 구성한다.
—
요소 방정식
각 요소에서
\[K^{(e)} U^{(e)} = F^{(e)} \tag{12.4.2-2}\]
강성행렬과 하중벡터는
\[K_{ij}^{(e)} = \int_{x_1}^{x_2} \left( a(x)\frac{d\phi_i}{dx}\frac{d\phi_j}{dx} + b(x)\phi_j\frac{d\phi_i}{dx} + \phi_i\phi_j \right) dx \tag{12.4.2-3}\]
\[F_i^{(e)} = \int_{x_1}^{x_2} f(x)\phi_i(x)\,dx\]
—
기준좌표계 적분 변환
좌표 변환:
\[\int_{x_1}^{x_2} f(x) dx = \int_{-1}^{1} f(x(\xi)) J d\xi \tag{12.4.2-4}\]
가우스 적분 적용:
\[\int_{-1}^{1} f(\xi) d\xi \approx \sum_{i=1}^{n_{GP}} w_i f(\xi_i)\]
여기서
\(\xi_i\) : Gauss-Legendre 적분점
\(w_i\) : 가중치
\(n_{GP}\) : 적분점 개수
—
도함수 좌표 변환
연쇄법칙:
\[\frac{d\phi}{dx} = \frac{d\phi}{d\xi} \frac{d\xi}{dx} = J^{-1} \frac{d\phi}{d\xi} \tag{12.4.2-6}\]
Jacobian:
\[J=\frac{dx}{d\xi} = \sum_{i=1}^{2} x_i \frac{d\phi_i}{d\xi} \tag{12.4.2-7}\]
—
최종 기준좌표계 요소 행렬
\[K_{ij}^{(e)} = \int_{-1}^{1} \left( a(x(\xi)) \frac{d\phi_i}{d\xi} \frac{d\phi_j}{d\xi} J^{-1} + b(x(\xi)) \phi_j \frac{d\phi_i}{d\xi} + \phi_i \phi_j \right) J \, d\xi \tag{12.4.2-8}\]
가우스 적분으로 계산하면
\[K_{ij}^{(e)} \approx \sum_{k=1}^{n_{GP}} w_k F_{ij}^{(e)}(\xi_k)\]
여기서
\[F_{ij}^{(e)}(\xi) = \left( a(x) \frac{d\phi_i}{d\xi} \frac{d\phi_j}{d\xi} J^{-1} + b(x) \phi_j \frac{d\phi_i}{d\xi} + \phi_i \phi_j \right) J\]
—
핵심 정리
실제좌표 적분 → 기준좌표 적분
Jacobian으로 미분 및 적분 변환
Gauss-Legendre 적분 적용
기준 요소에서 모든 계산 수행
전역 조립 후 선형 시스템 풀이
1. 미분방정식과 경계조건
다음 지배 방정식을 고려한다.
\[- \frac{d}{dx} \left( a(x)\frac{du(x)}{dx} \right) = f(x), \qquad x\in[0,1] \tag{12.4.3-1}\]
경계 조건:
\[u(0)=0, \qquad u(1)=g(1)\]
—
2. 형상함수와 도함수
기준 좌표 \(\xi \in [-1,1]\)에서
\[\phi_1(\xi)=\frac{1-\xi}{2}, \qquad \phi_2(\xi)=\frac{1+\xi}{2} \tag{12.4.3-2}\]
도함수:
\[\frac{d\phi_1}{d\xi}=-\frac{1}{2}, \qquad \frac{d\phi_2}{d\xi}=\frac{1}{2}\]
—
3. 물리좌표 변환
\[x(\xi) = x_1 \phi_1(\xi) + x_2 \phi_2(\xi) \tag{12.4.3-3}\]
정리하면
\[x(\xi) = \frac{x_2-x_1}{2}\xi + \frac{x_2+x_1}{2}\]
—
4. Jacobian
\[J = \frac{dx}{d\xi} = x_1 \frac{d\phi_1}{d\xi} + x_2 \frac{d\phi_2}{d\xi} = \frac{x_2-x_1}{2} \tag{12.4.3-4}\]
—
5. 물리좌표계 도함수 변환
연쇄법칙:
\[\frac{d\phi}{dx} = \frac{d\phi}{d\xi} \frac{d\xi}{dx} = J^{-1} \frac{d\phi}{d\xi} = \frac{2}{x_2-x_1} \frac{d\phi}{d\xi} \tag{12.4.3-5}\]
—
6. 강성행렬의 기준좌표 변환
요소 강성행렬:
\[K_{ij}^{(e)} = \int_{x_1}^{x_2} a(x) \frac{d\phi_i}{dx} \frac{d\phi_j}{dx} dx\]
기준좌표로 변환하면
\[K_{ij}^{(e)} = \int_{-1}^{1} a(x(\xi)) \frac{d\phi_i}{d\xi} \frac{d\phi_j}{d\xi} J^{-1} J \, d\xi \tag{12.4.3-6}\]
즉,
\[K_{ij}^{(e)} = \int_{-1}^{1} a(x(\xi)) \frac{2}{x_2-x_1} \frac{d\phi_i}{d\xi} \frac{d\phi_j}{d\xi} d\xi\]
—
7. 하중벡터 변환
\[F_i^{(e)} = \int_{x_1}^{x_2} f(x)\phi_i(x) dx\]
기준좌표 변환:
\[F_i^{(e)} = \int_{-1}^{1} f(x(\xi)) \phi_i(\xi) J \, d\xi \tag{12.4.3-7}\]
—
8. 가우스 적분을 이용한 근사
\[K_{ij}^{(e)} \approx \sum_{q=1}^{n_{GP}} w_q a(x(\xi_q)) \frac{2}{x_2-x_1} \frac{d\phi_i}{d\xi} \frac{d\phi_j}{d\xi} \tag{12.4.3-8}\]
\[F_i^{(e)} \approx \sum_{q=1}^{n_{GP}} w_q f(x(\xi_q)) \phi_i(\xi_q) \frac{x_2-x_1}{2}\]
—
핵심 정리
물리좌표 → 기준좌표 매핑
Jacobian으로 미분 및 적분 변환
강성행렬과 하중벡터를 기준좌표에서 계산
Gauss-Legendre 적분으로 수치 계산
1. 주어진 문제
\[- \frac{d}{dx} \left( a(x)\frac{du(x)}{dx} \right) = f(x), \qquad x \in [0,1] \tag{12.4.4-1}\]
경계 조건:
\[u(0)=0, \qquad u(1)=\cos(1)\]
계수 함수와 소스 항은
\[a(x)=e^{x}\]
\[f(x) = e^{x} \left( \cos(x) - 2\sin(x) - x\cos(x) - x\sin(x) \right) \tag{12.4.4-2}\]
정확해(Exact Solution):
\[u(x)=x\cos(x) \tag{12.4.4-3}\]
—
2. FEM 구현 절차
요소 개수 설정: num_elements
요소 길이 계산: \(h=1/\texttt{num\_elements}\)
요소 강성행렬 계산
요소 하중벡터 계산
전역 조립
경계조건 적용
선형 시스템 풀이
정확해와 비교
—
3. Python 코드
import numpy as np
import matplotlib.pyplot as plt
# ------------------------------
# Problem definition
# ------------------------------
def a(x):
return np.exp(x)
def f(x):
return np.exp(x)*(np.cos(x)
-2*np.sin(x)
-x*np.cos(x)
-x*np.sin(x))
def exact(x):
return x*np.cos(x)
# ------------------------------
# FEM parameters
# ------------------------------
num_elements = 10
num_nodes = num_elements + 1
x_nodes = np.linspace(0,1,num_nodes)
h = 1.0/num_elements
K = np.zeros((num_nodes,num_nodes))
F = np.zeros(num_nodes)
# ------------------------------
# Assembly
# ------------------------------
for e in range(num_elements):
x1 = x_nodes[e]
x2 = x_nodes[e+1]
# 2-point Gauss quadrature
xi = [-1/np.sqrt(3), 1/np.sqrt(3)]
w = [1,1]
for q in range(2):
# mapping
x = (x2-x1)/2*xi[q] + (x2+x1)/2
J = (x2-x1)/2
dphi = np.array([-1/2, 1/2])
phi = np.array([(1-xi[q])/2,
(1+xi[q])/2])
B = dphi * (2/(x2-x1))
Ke = a(x)*np.outer(B,B)*J*w[q]
Fe = f(x)*phi*J*w[q]
K[e:e+2,e:e+2] += Ke
F[e:e+2] += Fe
# ------------------------------
# Dirichlet BC
# ------------------------------
K[0,:]=0; K[0,0]=1
F[0]=0
K[-1,:]=0; K[-1,-1]=1
F[-1]=np.cos(1)
# ------------------------------
# Solve
# ------------------------------
U = np.linalg.solve(K,F)
# ------------------------------
# Plot
# ------------------------------
x_plot = np.linspace(0,1,200)
plt.plot(x_nodes,U,'o-',label='FEM')
plt.plot(x_plot,exact(x_plot),'r--',label='Exact')
plt.legend()
plt.show()—
4. 핵심 정리
기준좌표계에서 Gauss 적분 사용
Jacobian으로 좌표 변환
전역 강성행렬 조립
Dirichlet 경계조건 직접 대입
정확해와 비교 가능